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The two-dimensional antiferromagnetic S = 1/2 Heisenberg model with random bond dilution 
is studied using quantum Monte Carlo simulation at the percolation threshold (50% of the bonds 
removed). Finite-size scaling of the staggered structure factor averaged over the largest connected 
clusters of sites on L x L lattices shows that long-range order exists within the percolating fractal 
clusters in the thermodynamic limit. This implies that the order-disorder transition driven by bond- 
dilution occurs exactly at the percolation threshold and that the exponents are classical. This result 
should apply also to the site-diluted system. 
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During the past decade, questions related to the de- 
struction upon doping of the antiferromagnetic order 
in the high-T c cuprate materials have motivated ex- 
tensive studies of quantum critical phenomena in two- 
dimensional (2D) antiferromagnets 0. The 2D Heisen- 
berg model on a square lattice can be driven through 
an order-disorder transition by, e.g., introducing frus- 
trating interactions || or by dimerizing the lattice |Q. 
It has also been believed that a non-trivial phase tran- 
sition could be achieved by diluting the system, i.e., by 
randomly removing either sites |^-[| or bonds (nearest- 
neighbor interactions) p0| . The site-dilution problem is 
of direct relevance to the cuprates doped with nonmag- 
netic impurities, such as Zn or Mn substituted for Cu [[llj . 
Early numerical work indicated that the long-range order 
vanishes in the Heisenberg model with nearest-neighbor 
interactions when a fraction p* rj 0.35 of the sites are 
removed Q. Various analytical treatments have given 
results forp* ranging from 0.07 § to 0.30 §. These hole 
concentrations are below the classical percolation thresh- 
old p c \ w 0.407 Jl^], and hence the phase transition would 
be caused by quantum fluctuations. However, in a recent 
paper Kato et ai. reported quantum Monte Carlo simula- 
tions of larger lattices at lower temperatures than in pre- 
vious work and concluded that the critical site-dilution is 
exactly the percolation density; p* = p c \ Q . Based on 
their simulations, they also argued that the critical be- 
havior nevertheless is not classical, but that the fractal 
clusters at p c \ are quantum critical with algebraically de- 
caying correlation functions. This leads to non-classical 
critical exponents. Most surprisingly, the simulations in- 
dicated that the exponents depend on the magnitude of 
the spin. 

In this Letter we use a quantum Monte Carlo method 
to study the related bond-diluted system exactly at the 
percolation threshold, i.e., with 50% of the bonds ran- 
domly removed. In order to determine the nature of the 
ground state at this point — quantum critical, classically 
critical, or quantum disordered — we study the mag- 
netic properties of the largest clusters of connected sites 



on L x L lattices with L up to 18. We find clear evi- 
dence that these clusters, which in the thermodynamic 
limit are fractal with fractal dimension d = 91/48 [ pT[ , 
are antiferromagnetically ordered. This implies that the 
order-disorder transition driven by bond-dilution occurs 
exactly at the percolation threshold and that the critical 
exponents are classical. 

We have chosen to study the bond dilution problem 
because it is numerically more tractable than site dilu- 
tion — since the bond percolation threshold p c \ = 1/2 
it can be realized exactly on the finite lattices we work 
with. However, the fractal dimension and the critical 
exponents are the same for classical bond and site per- 
colation jlj|, and our conclusions should therefore hold 
true also for the site-diluted system. We argue that the 
reason for the disagreement with the previous results by 
Kato et al. [|l3) is that they did not use sufficiently low 
temperatures in their simulations — extremely low tem- 
peratures are required for reaching the ground state even 
for the relatively small system sizes we use here. 

We consider the 5=1/2 Heisenberg Hamiltonian on 
a square lattice with N = L x L sites; 

2N 

H = J2j(b)S i{b )-S j{b) . (1) 

6=1 

The bonds b connect nearest neighbor sites i(b),j(b) with 
interaction strength J(b) = J for N randomly selected 
bonds and Jib) — for the remaining N bonds. We 
use an efficient finite-temperature quantum Monte Carlo 
method based on the "stochastic series expansion" ap- 
proach [ ^6|p^ ] to study systems with L = 4, 6, . . . , 18. 
In order to reach the ground state, we successively in- 
crease the inverse temperature (3 = J/T until all quan- 
tities of interest have converged. We find that (3 as high 
as w 10 4 is required for the largest lattice we have con- 
sidered. Another important concern when studying ran- 
dom systems is the equilibration of the simulations. One 
would like to average over as many random configurations 
as possible within given computer resources. Ideally, one 
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would then perform only short simulations for each re- 
alization (typically, even a quite short simulation of a 
given configuration results in a statistical error smaller 
than the fluctuations between different configurations). 
However, there is a minimum length of a simulation set 
by the time needed to equilibrate it, and when carrying 
out short simulations it is important to have some way 
to verify that the correct equilibrium distribution indeed 
has been reached. We use the following scheme to check 
for both equilibration and temperature effects: For each 
bond configuration we carry out simulations at inverse 
temperatures (3 — 2 n L, n — 0,1, ... , n max . Starting with 
n = 0, we perform two runs for every (3, each with N e 
updating steps for equilibration and N m = 2iV e steps 
for measuring physical quantities (for the definition of a 
"step", see Ref. 0). The second run is a direct con- 
tinuation of the first one, so that the effective number of 
equilibration steps is four times that for the first run. An 
agreement between the results of these two runs is then 
a good indication that the simulation has equilibrated. 
For the subsequently lower temperatures (increasing n) 
we always start from the last Monte Carlo state gener- 
ated at the previous temperature. The convergence of 
the simulations using the /3-doubling procedure will be 
illustrated with som results below. 

In the thermodynamic limit, the system at p c \ will 
be spanned by infinite clusters of fractal dimension d = 
91/48 O]. The existence of a nontrivial (quantum) crit- 
ical point is determined by the magnetic properties of 
these clusters. If they are long-range ordered, the criti- 
cal point of the order-disorder transition driven by bond 
dilution will be exactly at p c \ = 1/2 as in the classical 
case, and the critical exponents will be the classical ones. 
If the fractal clusters are critical, i.e., their spin-spin cor- 
relations decay with a power-law as suggested by Kato 
et al. fl3| , the critical point is still at p c \ but the expo- 
nents are different and non-trivial. A quantum critical 
point with p* < p c \ would imply an exponential decay 
of the spin-spin correlations within the fractal clusters at 
Pel. In order to determine which of these scenarios ap- 
ply, we have calculated the magnetization squared of the 
largest connected cluster of sites for a large number of 
random bond configurations on lattices of different lin- 
ear size L. As L — > oo this procedure gives the ordered 
moment of the fractal clusters of interest. Denoting the 
largest cluster by C, the z-component of the staggered 
magnetization squared is given by 

-Hte- ir+ - sf ))' <2) 

where Nc is the number of spins in the cluster and the 
brackets indicate both the quantum mechanical expecta- 
tion value and an average over bond configurations. We 
also consider the full staggered structure factor 
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FIG. 1. Staggered magnetization squared for the largest 
cluster of connected spins on lattices with L = 18, averaged 
over 1000 bond configurations and normalized by the result 
obtained in the second simulation at the lowest temperature. 
The index n corresponds to inverse temperature according 
to f3 = 2 n L. The points for the second of the two equal-/3 
runs have been shifted to the right in order to enable a bet- 
ter comparison between the data sets. The inset shows the 
low-temperature data on a more detailed scale. 

^ = (^(E(- 1 )- 4 * fl ?) V ( 3 ) 

which involves all the spins of the lattice and was used 
by Kato et al. |f| to study the critical behavior of the 
site-diluted system. The number of random bond config- 
urations used in the averages presented here ranges from 
ps 10 3 for L = 18 to w 10 4 for L = 4. 

In Fig. [j] we show results illustrating the equilibration 
scheme. The disorder-averaged m c is graphed versus the 
index n (specifying the inverse temperature — 2 n L as 
described above) for L = 18. In order to reduce effects 
of fluctuations among bond realizations and more clearly 
show the relative effects of equilibration times and tem- 
perature, we have normalized the data to the result of 
the second run at the highest (3 used (n max = 9, cor- 
responding to j3 = 9216) and estimated the statistical 
errors using the bootstrap method fill . The number of 
equilibration and measurement steps were N c = 500 and 
N m = 1000. Within error bars, there are no differences 
between any of the equal-/? runs, and we therefore con- 
clude that the simulations are well equilibrated. As an 
additional check, for L = 16 and smaller we have also 
carried out simulations using N c — 250 and N c = 1000 
(keeping N m — 2N ). For N c = 250 we do see small but 
statistically significant differences between the equal-/? 
runs, but the results of the second run are always consis- 
tent with data obtained using the longer equilibrations. 
Hence, we believe that our results are free of detectable 
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FIG. 2. Disorder-averaged staggered structure factor vs in- 
verse temperature for systems of linear size L = 4, 6, . . . , 18 
(from top). The data has been normalized by the results ob- 
tained at the highest (3 for each system size {{3 = 2 9 L). 



effects of insufficient equilibration. All the results to be 
presented below were averaged using the second of the 
equal-/? runs only. 

Fig. |l| shows that very low temperatures are needed to 
converge to the ground state even for a lattice of rela- 
tively modest size. Since there are no statistically signif- 
icant differences between the L = 18 results at (3 = 4608 
and 9216, and the asymptotic approach to the ground 
state should be exponential, (3 — 9216 should give the 
ground state to within error bars. In Fig. || we show nor- 
malized results for S„ as a function of f3 for all the lattices 
we have studied. It is clear that using a fixed (3 = 500 
for all system sizes, as was done in Ref. [Q, leads to 
a significant systematic error (assuming that site-diluted 
systems are similarly affected by temperature, which can 
be expected up to some factor of order 1). Fig. shows 
that the relative deviation from the ground state grows 
rapidly with L at fixed (3 and is ss 3% for L = 18 at 
13 = 500. The largest lattice considered in Ref. |l3| was 
L = 48, for which our results suggest that (3 = 500 could 
lead to an error of more than 10%. The results to be 
discussed below are all for f3 = 2 9 L. 

The reason for the very high /3-values needed to con- 
verge to the ground state is most likely that localized 
moments can form in the irregular clusters of connected 
spins. These moments interact with each other with 
a strength which decreases rapidly with increasing sep- 
aration, thus leading to closely spaced energy levels. 
The typical level spacing should decrease faster with in- 
creasing L than the 1/L behavior suggested in Ref. 0. 
The temperature effects should be the largest exactly at 
the percolation threshold and for all hole concentrations 
they lead to an underestimation of the ordered moment. 
Hence, the result that there is a substantial order in the 
system close to the percolation threshold [|l3| should re- 




0.00 

0.00 0.05 0.10 0.15 0.20 0.25 



FIG. 3. Finite-size scaling of the disorder-averaged stag- 
gered magnetization squared of the largest connected clusters 
on L x L lattices. Statistical errors are smaller than the sym- 
bols. The curve is a fit to a cubic polynomial. 

main valid despite such effects. 

In order to definitely determine whether indeed p* = 
Pd and whether or not the exponents are the classical 
ones, we here study the lattice-size dependence of the 
cluster order parameter squared; Eq. (|^). As a finite- 
size scaling ansatz, we use a simple generalization of the 
known scaling law for the sublattice magnetization m for 
the pure 2D Hciscnbcrg model. In that case the lead- 
ing size-correction to m 2 is ~ 1/N 1 / 2 fll9| , which can be 
seen clearly in numerical data p0|| . Since the average 
number of spins in the fractal clusters of the diluted sys- 
tem depends asymptotically on L = N 1 / 2 according to 
(Nc) ~ L d , with d the fractal dimension quoted above, 
we here assume a leading size correction ~ 1/ ' L d l 2 . Fig. || 
shows our data for graphed versus this variable. The 
results appear to be consistent with the ansatz, although 
in order to fit all the data points we have to use a polyno- 
mial cubic in \jL d l 2 (a cubic polynomial is needed also 
to fit data for the pure 2D Heisenberg model, but the 
corrections to the linear behavior are much smaller |^0| ) . 
This fit has a \ 2 P er degree of freedom of w 0.6 and gives 
the full staggered magnetization Md = \Jim 2 c w 0.09 for 
the infinite fractal clusters (the factor 3 accounts for rota- 
tional averaging). Hence, the order on the ci-dimensional 
fractal lattice is as high as « 30% of the 2D staggered 
moment pp| ]. 

Finally, we discuss results for the staggered structure 
factor of the full lattice; Eq. (||). In Fig. |^ we graph 
In (5V) versus ln(L), which should show an asymptotic 
linear scaling behavior. There is a clear upward curvature 
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as the lattice size increases, and it is clear that we are still 
far from the scaling regime. A curvature is also seen in 
the data presented by Kato et al. for small systems [[l3| , 
but for larger sizes a linear behavior was discerned. As we 
have discussed above, the structure factor calculated for 
large lattices in Ref. p| was likely substantially under- 
estimated due to temperature effects, and the apparent 
non-classical scaling behavior is then an artifact. Since 
the results in Fig. || show that there is long-range order in 
the fractal clusters, the growth of must asymptotically 
be given by classical percolation theory, i.e., S n ~ L 2d ~ 2 , 
where 2d — 2 = 43/24 Hi. The very large corrections 
to scaling evident in Fig7g| can be understood as result- 
ing from a significant reduction with increasing L of the 
staggered structure factor per site of the fractal clusters 
(m 2 -,), as seen in Fig. 0. Classically, is independent 
of system size and the scaling of S„ is therefore deter- 
mined solely by the increase in size of the clusters as L 
increases. In the quantum case, the size dependence of 
SV is dominated by this geometric effect only for sys- 
tem sizes sufficiently large for the relative size-correction 
to m 2 -, to be small. For our largest lattice, the relative 
size correction is still more than a factor four. Consid- 
ering the extremely low temperatures required to study 
the ground state of large lattices, it will be very difficult 
to numerically observe the asymptotic classical scaling 
regime for S = 1/2. 

In summary, we have presented numerical results show- 
ing that the fractal S — 1/2 Heisenberg clusters at the 
classical bond-percolation density have long-range anti- 
ferromagnetic order. This implies that the order-disorder 
transition driven by bond-dilution occurs exactly at the 
percolation density and that the critical exponents arc 
classical. This should hold true for any spin S, and 
since classical bond and site percolation are equivalent 
in terms of the fractal dimension and the critical expo- 
nents jl5|) , our conclusions should apply also to the site- 
diluted model. The quantum mechanical corrections to 
the asymptotic scaling behavior are very large for lattice 
sizes that can be studied with today's computers, mak- 
ing direct observation of the critical behavior difficult. 
We have also discussed temperature effects and shown 
that extremely low temperatures are needed to study the 
ground state, likely due to the presence of weakly inter- 
acting localized moments. 

It would be interesting to study a diluted system in- 
cluding frustration. If the clean system is already close to 
a quantum critical point, non-magnetic impurities may 
be able to drive it into a quantum disordered phase. 
The fact that Zn or Mn doping of the cuprates destroys 
the long-range order well before the classical percolation 
threshold [JllJ clearly indicates that these materials can- 
not be described by a randomly diluted Heisenberg model 
with nearest-neighbor interactions only. 

This work was supported by the NSF under grant 
No. DMR-97-12765. Some of the calculations were car- 
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FIG. 4. The logarithm of the staggered structure factor vs 
the logarithm of the system size. The line has the slope 43/24 
expected in the asymptotic classical scaling regime. 



ried out using the Origin2000 system at the NCSA. 
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